!-------------------------------------------------!
!               MODULE HAMILTONIAN		  !
!-------------------------------------------------!
! This module produces the hamiltonian needed     !
! in the program		                  !
!-------------------------------------------------!

  module hamiltonian
  
  use constants
  use control 
  use matrixmod
    
  implicit none
    
  private 
    
!-------------------------------------------------!
! Variables declaration	          		  !					
!-------------------------------------------------!
!-------------------------------------------------!
! Public variables  				  !	
!-------------------------------------------------!
    
  real(dp), public, allocatable :: H(:,:), Solap(:,:)
  real(dp), public, allocatable :: r_atoms(:,:), self_orbital(:)
  integer, public :: N_atoms, N_self, N_deco
  integer, public, allocatable :: atom_number(:), deco_orbital(:)
  character(len=2),public, allocatable :: atom_type(:) 

!-------------------------------------------------!
! Public subroutines	 	        	  !					
!-------------------------------------------------!
    
  public :: h_read,  S_read, write_H_matrix, coords_read, self_energy_read_from_axis, self_energy_read_from_file
  public :: decoherence_read

  contains

!-------------------------------------------------!
! 		H_READ 	 		  	  !					
!-------------------------------------------------!
! This reads the hamiltonian from a file          !
! Hamiltonian must always be in (eV) 	          !
!-------------------------------------------------!

  subroutine h_read()
    
    implicit none
    integer :: i,j,cont,Naux
    real(dp) :: prom,conv_units


    if(Banner_logic.EQ.'S') write(*,*)'readind hamiltonian ...'
   
    open(1,file='hamiltonian.dat')
    read(1,*)Naux
    write(*,*)N,Naux  
    if(Naux.NE.N) stop'Error: Matching error between coordinates and hamiltonian' 

    allocate(H(N,N))
      
    H=0.0_dp
      
    prom = 0.0_dp
      
    cont = 0 
    do i=1,N
      do j=1,N 
        read(1,*)H(i,j)
        cont = cont + 1        
      enddo
      prom = prom + H(i,i)
    enddo

    if(Banner_logic.EQ.'S') write(*,*)'DimensionOfHamiltonianMatrix NxN',cont
    prom = prom/real(N)

    if(method_logic.EQ.'HUCK')conv_units = 1.0_dp
    if(method_logic.EQ.'DFTB')conv_units = hartreetoev
    if(method_logic.EQ.'STCH')conv_units = 1.0_dp
    if(method_logic.EQ.'SIES')conv_units = hartreetoev/2.0_dp !Ritbergs

    do i=1,N
      do j=1,N
        if(j.GT.i)H(i,j)=H(j,i)
  	H(i,j)=conv_units*H(i,j)
      enddo
    enddo
	
    close(1)
      
  end subroutine h_read

!-------------------------------------------------!
! 		Coords_read                     		  !					
!-------------------------------------------------!
! This reads the coordinates and the atom type	  !
!-------------------------------------------------!
  subroutine coords_read()
    
    implicit none
    integer :: i,j
    integer :: N_orb

	    
    open(1,file='coords.xyz')		
    if(Banner_logic.EQ.'S') write(*,*)'reading coordinates ...'
    read(1,*)N_atoms
    read(1,*)
    if(Banner_logic.EQ.'S') write(*,*)'NumberOfAtoms',N_atoms
    allocate(atom_type(N_atoms))
    allocate(r_atoms(N_atoms,3))
    allocate(atom_number(N_atoms))
      
    r_atoms=0.0_dp
      
    do i=1,N_atoms
      read(1,*)atom_type(i),r_atoms(i,1),r_atoms(i,2),r_atoms(i,3)
    enddo
  
    do i=1,N_atoms
      do j=1,NZ
        if(atom_type(i).EQ.NAME(j))atom_number(i)=j
      enddo
    enddo
                 
    close(1)
      
    N_orb=0       
      do i=1,N_atoms
        N_orb=N_orb + NORB_METHOD(atom_number(i))
      enddo
    
    if(Banner_logic.EQ.'S') write(*,*)'NumberOfOrbitals =',N_orb

    N=N_orb
      
  end subroutine coords_read
    
!-------------------------------------------------!
! 		Self_energys_positions            		  !					
!-------------------------------------------------!
! This reeds the position of the Selfs_energys	  !
!-------------------------------------------------!
  subroutine self_energy_read_from_axis()
   
  implicit none
      
    integer :: i,j,k,cont
    integer, allocatable :: self_atom(:,:)
		character(20), allocatable :: self_energy_pos(:)
  
    N_self=N_atoms 
    allocate(self_atom(N_self,2))
    allocate(self_orbital(N))
    allocate(self_energy_pos(N_atoms))

      
    self_orbital = 0
    self_atom = 0
      
    if(Banner_logic.EQ.'S') write(*,*)'reading selfs energys positions ...'
     
    cont=0
    self_energy_pos = '1 1 1'
  
    do i=1,N_atoms
      if(Conducting_axis.EQ.'x')then
        if(r_atoms(i,1).lt.LSEP)then
          cont=cont+1
          self_atom(cont,1)=i
          self_atom(cont,2)=-1
        endif
        if(r_atoms(i,1).gt.RSEP)then
          cont=cont+1
          self_atom(cont,1)=i
          self_atom(cont,2)=1
        endif 
      endif  
      
      if(Conducting_axis.EQ.'y')then
        if(r_atoms(i,2).lt.LSEP)then
          cont=cont+1
          self_atom(cont,1)=i
          self_atom(cont,2)=-1
        endif
        if(r_atoms(i,2).gt.RSEP)then
          cont=cont+1
          self_atom(cont,1)=i
          self_atom(cont,2)=1
        endif
      endif       
      
      if(Conducting_axis.EQ.'z')then
        if(r_atoms(i,3).lt.LSEP)then
          cont=cont+1
          self_atom(cont,1)=i
          self_atom(cont,2)=-1
          self_energy_pos(i)='0 0 0'
        endif
        if(r_atoms(i,3).gt.RSEP)then
          cont=cont+1
          self_atom(cont,1)=i
          self_atom(cont,2)=1
          self_energy_pos(i)='0 0 0'
        endif
      endif  
    enddo
      
    N_self = cont
     
    if(Banner_logic.EQ.'S') write(*,*)'NumberOfSelfsEnergys',N_self
    
    !Esto es un simple algoritmo de llenado de orbitales
      
    do i=1,N_self
      cont=0
      do j=1,self_atom(i,1)-1
        do k=1,NORB_METHOD(atom_number(j))
          cont=cont+1
        enddo
      enddo
      do k=1,NORB_METHOD(atom_number(self_atom(i,1)))
        cont=cont+1
        self_orbital(cont)=self_atom(i,2)
      enddo 
    enddo


    write(999,*)N_atoms
		do i=1,N_atoms
			write(999,*)i,self_energy_pos(i)
		enddo


  end subroutine self_energy_read_from_axis

!-------------------------------------------------!
!  Read the self energy from file 		  !
!-------------------------------------------------!

   subroutine self_energy_read_from_file()
   
      implicit none
      
      integer :: i,j,k,cont
      integer, allocatable :: self_atom(:)
      real, allocatable :: loads(:,:)
      
      open(1,file='self_energy.inp',status='OLD')
      read(1,*)N_self
          
      allocate(self_atom(N_self))       !four columns: i , ns, np, nd
      allocate(self_orbital(N))
      allocate(loads(N_self,3))
      
      self_orbital = 0
      self_atom = 0
      
      if(Banner_logic.EQ.'S') write(*,*)'reading selfs energys positions ...'
      
      do i=1,N_self
       read(1,*)self_atom(i),loads(i,1) ,loads(i,2),loads(i,3)
       write(*,*)self_atom(i),loads(i,1) ,loads(i,2),loads(i,3)
      enddo
      close(1)
 
     if(Banner_logic.EQ.'S') write(*,*)'Atoms with SelfsEnergys',N_self
      

     if(METHOD_logic.EQ.'STCH'.or.METHOD_logic.EQ.'HUCK'.or.METHOD_logic.EQ.'DFTB')then
       do i=1,N_self
        cont=0
        do j=1,self_atom(i)-1
         do k=1,NORB_METHOD(atom_number(j))
          cont=cont+1
         enddo
        enddo


        do k=1,NORB_METHOD(atom_number(self_atom(i)))
          cont=cont+1
          if(NORB_METHOD(atom_number(self_atom(i))).EQ.1) self_orbital(cont)=loads(i,1) !s

          if(NORB_METHOD(atom_number(self_atom(i))).EQ.4) then
            if(k.EQ.1) self_orbital(cont)=loads(i,1) !s
            if(k.GT.1) self_orbital(cont)=loads(i,2) !p
          endif

          if(NORB_METHOD(atom_number(self_atom(i))).EQ.9) then
            if(k.EQ.1) self_orbital(cont)=loads(i,1) !s
            if(k.GT.1.and.k.LT.5) self_orbital(cont)=loads(i,2) !p
            if(k.GE.5) self_orbital(cont)=loads(i,3) !d
          endif	
        enddo 
       enddo

     endif


     if(METHOD_logic.EQ.'SIES')then
       do i=1,N_self
         cont=0
         do j=1,self_atom(i)-1
           do k=1,NORB_METHOD_s(atom_number(j))
             cont=cont+1
           enddo
           do k=1,NORB_METHOD_p(atom_number(j))
             cont=cont+1
           enddo
           do k=1,NORB_METHOD_d(atom_number(j))
             cont=cont+1
           enddo
         enddo

         do k=1,NORB_METHOD_s(atom_number(self_atom(i)))
	   cont=cont+1
           self_orbital(cont)=loads(i,1) !s
!	   if(k.EQ.1)self_orbital(cont)=0.2 !para poder seleccionar el primer orbital
         enddo	
        
         do k=1,NORB_METHOD_p(atom_number(self_atom(i)))
   	    cont=cont+1
            self_orbital(cont)=loads(i,2) !p
         enddo	

         do k=1,NORB_METHOD_d(atom_number(self_atom(i)))
 	   cont=cont+1
           self_orbital(cont)=loads(i,3) !d
         enddo
       enddo	
     endif   

    
    end subroutine self_energy_read_from_file


!-------------------------------------------------!
! 		Decoherence Read - betha 	  !			
!-------------------------------------------------!
! This reads the decoherence from file            !
!-------------------------------------------------!

   subroutine decoherence_read()
   
      implicit none
      
      integer :: i,j,k,cont
      integer, allocatable :: deco_atom(:,:)

      open(1,file='decoherence.inp')
      read(1,*)N_deco
          
      allocate(deco_atom(N_deco,2))
      allocate(deco_orbital(N))
      
      deco_orbital = 0
      deco_atom = 0
      
      if(Banner_logic.EQ.'S') write(*,*)'reading decoherece positions ...'
      do i=1,N_deco
       read(1,*)deco_atom(i,1),deco_atom(i,2)
!       write(*,*)self_atom(i,1),self_atom(i,2)
      enddo
      close(1)
    
     if(Banner_logic.EQ.'S') write(*,*)'NumberOfDecoherencePoints',N_deco
      
     do i=1,N_deco
      cont=0
      do j=1,deco_atom(i,1)-1
       do k=1,NORB_METHOD(atom_number(j))
        cont=cont+1
       enddo
      enddo
      do k=1,NORB_METHOD(atom_number(deco_atom(i,1)))
        cont=cont+1
        deco_orbital(cont)=deco_atom(i,2)
      enddo 
     enddo

    end subroutine decoherence_read

!-------------------------------------------------!
! 		H_WRITE		                  !					
!-------------------------------------------------!
! This writes the hamiltonian matrix  	          !
!-------------------------------------------------!
  
  subroutine write_H_matrix()
   
    implicit none
    complex(dp), allocatable :: Aux(:,:)
    
    allocate(Aux(N,N))
    Aux=H*cmplx(1.0_dp,0.0)
    call write_matrix(Aux)      
 
  end subroutine write_H_matrix

!-------------------------------------------------!
! 		S_READ		 		  !					
!-------------------------------------------------!
! This reads the overlap matrix from a file       !
!-------------------------------------------------!

  subroutine S_read()
    
    implicit none
    integer :: i,j,Naux, cont
    complex(dp), allocatable  :: Solap1(:,:)
      	
    if(Banner_logic.EQ.'S') write(*,*)'reading overlap matrix ...'	
    open(1,file='overlap.dat')		
    read(1,*)Naux
    if(Naux.NE.N) stop'The dimension overlap matrix does not match the dimension of the hamiltomian matrix'
  
    allocate(Solap(N,N))
    ALLOCATE(Solap1(N,N))       
      
    Solap = 0.0_dp
  
    do i=1,N
      do j=1,N
        read(1,*)Solap(i,j)
      enddo
    enddo
      
    cont=0
 	
    do i=1,N
      do j=1,N
	if(j.GT.i)Solap(i,j)=Solap(j,i)
	  cont=cont+1
      enddo
    enddo
       
    if(Banner_logic.EQ.'S') write(*,*)'DimensionOfOverlapMatrix=',cont 
       
    Solap1=Solap

    close(1)
      
  end subroutine S_read

!-------------------------------------------------!
! End of the module 				  !
!-------------------------------------------------!

 end
